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Abstract 

We present a novel method for the calculation of the fractal dimension of 
boundaries in dynamical systems, which is in many cases many orders of 
magnitude more efficient than the uncertainty method. We call it the Out- 
put Function Evaluation (OFE) method. The OFE method is based on an 
efficient scheme for computing output functions, such as the escape time, on 
a one-dimensional portion of the phase space. We show analytically that 
the OFE method is much more efficient than the uncertainty method for 
boundaries with D < 0.5, where D is the dimension of the intersection of the 
boundary with a one-dimensional manifold. We apply the OFE method to 
a scattering system, and compare it to the uncertainty method. We use the 
OFE method to study the behavior of the fractal dimension as the system's 

dynamics undergoes a topological transition. 
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Motivation. - In dynamical systems with two or more well-defined asymptotic states 
(examples are found in systems with attracting or repelling sets), the boundary in phase 
space separating initial conditions corresponding to distinct final states (that is, belonging 
to different basins) may have a fractal structure; in this case one says that the system 
has a fractal basin boundary This means that sets of points in phase space which 
undergo very different time evolutions are mixed in a complex way in all scales. Fractal 
basin boundaries are found in many important physical systems, such as in astrophysics || , 
scattering systems |§, systems with escapes [|J, dissipative systems ||, etc. The occurrence 
of fractal basin boundaries implies a great sensitivity of the long-time evolution of the 
system to perturbations in the initial conditions. This sensitivity to initial conditions is 
characterized by the box-counting dimension d of the basin boundary, which is interpreted 
as a measure of the degree of uncertainty about the final fate of a system with fractal basin 
boundaries 0. The box-counting dimension is one of the most important quantities for 
characterizing these boundaries. Denoting the dimension of the total phase space by d ps , d 
satisfies in general d ps — l<d< d ps . d = d ps — 1 for regular boundaries, and d > d ps — 1 for 
fractal boundaries. 

Because of its fundamental physical significance, it is very important to have efficient 
numerical methods for calculating d with good precision. The best method known so far is 
the uncertainty method, which is based on a direct exploitation of the final state uncertainty 
reflected in d || . The uncertainty method is very efficient for high values of d (that is, for d 
close to d ps ), but it is inefficient for low values of d (close to d ps — 1), when the basin boundary 
departs only a little from a smooth manifold. This is because the number of initial points 
whose orbits are integrated for a calculation of d in this method for a one-dimensional set of 
initial conditions behaves as e D_1 , where D is the reduced dimension D = d — d ps + 1, and e 
is the smallest scale used in the computation, which defines the precision in the calculation: 
the smaller e is, the higher the precision will be. D is the dimension of the intersection of 
the boundary with a generic one- dimensional segment in phase space, and is bounded by 
< D < 1. We see that for D close to 0, this method is very inefficient, because in this case 
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the number of integrated orbits grows rapidly as e decreases. Motivated by this, we introduce 
in this Letter a new method for the numerical calculation of the box-counting dimension 
that is highly efficient for small values of D, and it is in this respect complementary to the 
uncertainty method. We call our method the Output Function Evaluation (OFE) method. 
Specifically, we show that the number of integrated orbits in the OFE method scales as <r D . 
We then apply the method to a scattering system and compare its results and performance 
with the uncertainty method. We show that the two methods give the same result (as of 
course they should), but that for low D our method is several orders of magnitude more 
efficient. Finally, we apply the OFE method to this scattering system and we show that the 
fractal dimension shows a characteristic behavior at a critical energy for which the invariant 
set suffers a topological transition. 

Method. — We start with a brief exposition of the uncertainty method. Consider a one- 
dimensional set C in the phase space. Take now a pair of points in C separated by a distance 
/, in a random position in C. The probability P(l) that the two points belong to different 
basins scales as P(l) ~ l x ~ D . The uncertainty method amounts to a direct calculation of 
P(l) for many values of I, thereby obtaining D. This is done by choosing randomly many 
pairs of /-separated points, and numerically integrating their corresponding trajectories. For 
a large number of pairs, the fraction of pairs which evolve to distinct final states for a given 
separation / should approach P(l). Doing this for several values of /, one can find D by 
fitting P{1) to a power law. The precision of the resulting D obtained in this way depends 
ultimately on the smallest value for / used in the computation, which we denote by e. For e 
small enough, the total number of initial conditions N unc integrated in the computation of 
the uncertainty method scales as 1/P(e), that is, 

N unc (e) ~ e D ~\ (1) 

Since the computation time is roughly proportional to N unc , one sees that although being 
very efficient for D close to one, this method is definitely inefficient for D near zero. 

The reason why the uncertainty method is inefficient for D close to zero is that for low 
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D the volume occupied by the boundary for a given resolution e (with e small) is very small, 
and the vast majority of pairs of initial conditions will fall on the same basin. Since one 
needs to find a minimum number of pairs of points belonging to different basins in order 
to have a reasonable statistics for P, this causes N unc to become very large, and to grow 
very rapidly as e decreases. To calculate low D values, we need a method that "focuses its 
attention" on the vicinity of the boundary, and concentrates its evaluations there. This is 
exactly what the OFE method does, and we explain it next. 

Our method is based on the computation of suitable output functions of the system; these 
are functions relating the values of particular variables of the system after it reaches one of 
the final states (e.g., after it converges to a neighborhood of an attractor in a dissipative 
system) to the initial conditions that led to that time evolution. One example is the time r 
it takes to reach the final state; another example is the deflection angle 0(6) as a function 
of the impact parameter b in scattering systems. We restrict ourselves to output functions 
whose domain is a one-dimensional sub-manifold of the phase-space which intersects the 
basin boundary. If the boundary is fractal, so is its intersection with this sub-manifold. 
The output functions mirror the geometrical structure of the basin boundary, and if the 
boundary is fractal, so are the output functions, and in this case they have a fractal set of 
singularities with the same dimension D as the basin boundary. 

Although our method works with many choices of output functions, in order to clarify 
the exposition, from now on we assume for definiteness that we are dealing with a scattering 
system, and we choose the output function to be the deflection angle 0(6). We have to 
calculate with a resolution given by e, with e being very small. The straightforward 
method of laying on a given interval of b an e-size grid and calculating for the points 
on the grid is not good, since the number of integrations goes as e _1 , which is even more 
inefficient than the uncertainty method. To improve this, we use a variable-sized grid. The 
idea is to adjust the size of the grid on b (the stepsize) so that the oscillations in are well 
resolved, with the minimum size of the grid given by the resolution e. For values of b away 
from the boundary, is smooth, and the grid size can be large, whereas for b close to the 

4 



boundary, is very steep and typically shows very wild oscillations; in this latter case, the 
grid size needs to be small to resolve 0. The steepness of for a given b is measured by 
the modulus of its derivative \d(p/db\, and we choose the gridsize A to be proportional to 
l/\d(j)/db\, with the constraint A > e. In this way, most of the computing will happen near 
the fractal region of 0, that is, near the basin boundary. 

The method is implemented as follows. For a given b interval (bi n ,bfi n ), is calculated 
sequentially for a set of b values 60, &2, • • • , bj, ■ ■ •, with bj < bj + \. We proceed by first 
choosing 6 = hn, b\ = b. in + e, and by integrating initial conditions corresponding to 60 and 
b\, we compute 0(&o) and 0(&i), which we denote by 0o and 0i, respectively. Generally, we 
use the notation <pj = <p(bj). Now from b and &i we obtain b 2 from 6 2 = &i + Ai, where the 
stepsize Aj is given by 

C75 if c ^ ^ A ma x 
Aj= e, if^<e (2) 

^■maxi if Cj ^max 

with 

where A max , 5, and a are constant parameters. d<f)(bj)/db is the derivative of calculated 
at b = bj. The idea is that the stepsize Aj be chosen so that J+1 — <f>j w 5, to a first- 
order approximation. In other words, the stepsize is chosen so that the variation of <fi from 
one point to the next is kept approximately constant; this is the key idea of our method. 
However, we do not allow the stepsize to grow too much from one step to the next: Eq. ([3D 
ensures that Aj/Aj-i < a, with a > 1 giving the constraint on the growth of the stepsize. 
This avoids problems near extrema, where d<f)(b)/db = and the first-order estimate of 
4>j + i — 4>j is not valid. Also, Aj is restricted to be within the interval (e, A max ). We use the 
two-point approximation for the derivative d<fi(bj) / db: 



Now from b%, we calculate 02, and obtain 63 through 63 = 62 + A2, and so on. The compu- 
tation is stopped when we reach step N for which > bfi n . 

Once we have calculated <j)(b) by this procedure, the fractal set of singularities is given (to 
resolution e) by points bj such that \bj— bj-\\ > f3, with (3 being a parameter satisfying (3 > S. 
This means that we pick the points bj corresponding to regions where the output function 
shows oscillations on scales smaller than e. The results of the method are independent of 
the value chosen for (3. 

Now that we are in the possession of a set of points M approximating the boundary to 
resolution e, we proceed to calculate the fractal dimension. We could use a direct implemen- 
tation of the definition of the box-counting dimension, but we prefer to use a more powerful 
method, described in [7||§, which we explain briefly now. For each point 6j in M, we count 
the number rii(l) of points in M that lie within a distance / of 6j. It can then be shown |7j 
that the average of l/nj(/) over all points in M scales with I as 

where D is the reduced dimension. From Eq. (^j), we obtain D by calculating (l/n(/)) for 
many values of I (with I small), and fitting the results to a power law. 

We now estimate how the number of integrations Nofe of the OFE method scales with 
the resolution e. By construction, in the OFE method most of the integrations are performed 
for points in the vicinity of singularities in the output function, where the function is the 
steepest. For a small enough e, N OFE is therefore proportional to the number of singularities 
that are resolved with resolution e; but by the definition of the box-counting dimension, this 
number is proportional to e~ D . Therefore, we have the result 

Nofe ~ e~ D '. (6) 

For D close to zero, Nqfe grows very slowly, and the OFE method is much more efficient 
than the uncertainty method; the opposite holds for D close to one. To better compare the 
two methods, we define /(e) to be the ratio of Nqfe and N unc . From Eqs. (|I|) and (0), we 
have: 
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/(e) = 2pi ~ . (7) 
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We see that / — > for e^0ifD<0.5. This means that for D < 0.5 (and e small enough) 
the OFE method is more efficient than the uncertainty method, and it becomes ever more 
so as e decreases. In fact, we will see in the example that follows that the difference in 
efficiency can be of many orders of magnitude. On the other hand, if D > 0.5, / — > oo for 
e — > 0; in this case, the uncertainty method is the more efficient one. 

Example. - - We exemplify our method with a Hamiltonian scattering system with two 
degrees of freedom, described by a potential function V(x,y), where V is required to be 
highly localized around the origin. To exemplify our results, we use a potential that is a 
superposition of three repulsive Gaussian hills: 
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y(x,y) = 5><exp(-r?/2o?), (8) 
where 

rj = (x-x i ) 2 + (y-y i ) 2 . (9) 

(x{, yi) are the coordinates of the centers of the three hills. Vi and Oi are constants, and give 
the height and the spread of each hill, respectively. Potentials of the form (^j) are paradigms 
of chaotic scattering, and have been extensively studied [IT|. For our example, we choose 



the parameters to be: X\ = —x 2 — 4, yi — y 2 — 0, x 3 = 0, y 3 = 2; V\ — V% — 10, V 3 = 1; 
o\ = cr 2 = 0.75, cr 3 = 0.325. For this potential, there is a transition from regular to chaotic 
scattering as the energy of the incoming particle drops below a critical energy E c , with 
E C >1. 

We now proceed to apply the OFE method to this system. For the output function we 
choose the deflection function <p(b), calculated for a one-dimensional set of initial conditions 
on the segment y = —10, < x < 4, with initial velocity parallel to the y axis. The initial 
position is sufficiently far away from the origin so that the particle can be considered to be 
initially in free motion, and the velocity is fixed by the energy constraint Vq = \ / 2E. In this 



case, the impact parameter b is simply the x coordinate. Each initial condition is integrated 
until its distance from the origin becomes greater than 10, when it can be considered to 
be in free motion again. The deflection angle is found by integrating the quadrature 
9 = (xv y — yv x )/y / x 2 + y' 2 along the trajectory, with being given by the value of 9 after 
the particle is scattered. 

We now calculate 0(6) by the OFE method, with e = 10~ 10 , a = 2, 5 = 0.03, and 
A max = 10~ 3 , for E = 0.95. We then find the approximation M to the fractal set of 
singularities of as explained above, using (3 = 0.5. Next we calculate the fractal dimension 
from M, using the method described above []7|||. The result is shown in Fig. la, which is a 
plot of the logarithm of (l/n(l)) as a function of I. The points clearly define a straight line, 
and from its inclination we get the dimension: D = 0.238 ± 0.002. We have also calculated 
D by the uncertainty method. In Fig. lb, the fraction g(l) of /-separated pairs of initial 
conditions with differing by at least (3 = 0.5 are plotted against I, in a log-log plot. For 
each I, we keep integrating pairs of points until 100 pairs have been found with differing 
by (3 or more. D is found from a least-square fit of the points in Fig. lb and using Eq. 
(0) to be D = 0.24 ± 0.02. The results of the two methods agree, as they should, but the 
result of the OFE method is ten times more accurate than that of the uncertainty method, 
even though the number of integrated points in the OFE calculation was only about 26000, 
compared to about 1.5 x 10 6 integrations that were necessary in the uncertainty method. 
Note that in Fig. 2b the function g was not extended to e = 10~ 10 , because that would 
require an unreasonable number of integrations. From Eq. ([!]) we can estimate the number 
of integrations N unc necessary to extend g down to I = e; using the value 0.24 for D, we 
find N unc to be over 10 10 , which is prohibitively large, and is six orders of magnitude larger 
than the number of integrations in the OFE method. This shows the superiority of the OFE 
method over the uncertainty method for D < 0.5. 

The potential (§) undergoes a topological change in the phase space as the energy drops 
below E = 1, due to the appearance of a new forbidden region around hill 3. At this energy 
there is a sudden change (a metamorphosis) in the topological structure of the invariant set 
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IIJ. We expect this metamorphosis to imply a characteristic behavior of D in the vicinity 



of E — 1. To test this, we apply the OFE method to obtain D as a function of E for E 
close to 1. The result is shown in Fig. 2. We see that D has a minimum at E — 1, and 
D(E) exhibits a cusp at this energy. Notice that a high accuracy in the computation was 
necessary to resolve the behavior of D at the metamorphosis. The calculation of Fig. 3 
using the uncertainty method would require a prohibitively large amount of computation to 
come up with the same results. 

We note that in the particular case of scattering in two dimensions, it can be shown that 
all the thermodynamical quantities associated with the fractal invariant set can be obtained 
from the time delay function [HJ. This means that we can use the OFE method to compute 
any such quantity, including for example the topological entropy. 

In summary, we have presented the OFE method for the calculation of the fractal dimen- 
sion which is much more efficient than the uncertainty method for D < 0.5. We illustrated 
the method with a scattering system, and we have shown for that case that our method is 
many orders of magnitude more efficient than the uncertainty method. We used the OFE 
method to show that the fractal dimension displays a characteristic behavior at a topological 
transition of the well-known three-hill system, which would have been very difficult to be 
resolved by the uncertainty method. 

This work was sponsored by FAPESP and ONR(Physics). 
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FIG. 2. Reduced dimension D as a function of the energy E. 
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